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ABSTRACT 

A possible inconsistency arising when a radiative cooling term is incorporated in a 
finite resolution self-gravitating hydrodynamic code is discussed. The inconsistency 
appears when the heating-cooling balance within the cooling and collapsing gas cloud 
is broken near the resolution limit of a numerical code. As the result, the cooling 
time of a fluid element increases enormously, leading to the unphysical conclusion 
that the fluid element does not cool and is therefore stable against the collapse. A 
special cooling consistency condition is introduced which approximately restores the 
heating-cooling balance and leads to a numerical solution that closely mimics the exact 
(infinite resolution) solution. 
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1 INTRODUCTION 



Cosmological hydrodynamic simulations are now playing 
larger and larger role in developing the theory of galaxy for- 
mation. They are successfully incorporating more and more 
complex physical processes that are important in the process 
of galaxy formation, such as dark matter evolution, gas dy- 
namics, radiative processes, star formation etc. Simulations, 
therefore, are challenged to solve a very complex system of 
nonlinear equations incorporating physical processes with 
widely different characteristic time-scales. However, one usu- 
ally does not solve the whole system of equations by a sin- 
gle numerical technique (simply because there are no such 
methods invented yet), but by combining various numeri- 
cal techniques for separately solving for dark matter, gas 
dynamics, radiative cooling etc. The advantage of this ap- 
proach is that one can use known numerical methods that 
have been developed to separately simulate each of physical 
processes incorporated into a larger simulation. However, 
the pay off for using different methods for solving parts of 
the whole problem is that those methods are not necessarily 
consist ent with each other; in other words, straightforward 
combin 



dynamics ought to be satisfied in order to produce physically 
sensible solutions. 

However, the requirement of consistency between nu- 
merical methods solving separate parts of the whole system 
of equations is a general one, and applies not only to the 
case of gravity and gas dynamical solvers. In particular, in 
this paper I demonstrate how a cooling consistency condi- 
tion can emerge when radiative cooling is incorporated into 
a self-gravitating hydrodynamic code. 

Radiative cooling is the crucial mechanism responsible 
for condensation and collapse of baryonic gas into galaxies, 
and an adequate treatment of it in numerical simulations is 
required in any realistic cosmological gas dynamical simula- 
tion except, perhaps, in a simulation of clusters of galax- 
ies. Yet, since radiative cooling can often introduce into 
the solution the time-scales that are orders of magnitude 
shorter than the characteristic time-scale for the gas evo- 
lution (the sound crossing time), no numerical scheme has 
been invented to solve the whole system of equations includ- 
ing gas dynamics and cooling together, but rather two differ- 
ent numerical methods, one for gas dynamics, and another 
for cooling, are combined together t o incorporate radiative 
cooling into the gas dynamical code (Cen 1992; [Katz, Hem 1 



ferent sub-sets of the whole system of nonlinear equations 
may not lead to a techniques that solves the whole system. 
An example of how direct merging of a gas dynamical solver 
with a gravity solver can lead to a code that produces un- 
physical results (energy i s not balanced locally) is given in 
Gnedin & Bertschinger (1996). Therefore, a special gravi- 
tational consistency condition between the gravity and gas 



ng of two of numerical techniques that solve two dif- quist, fc Weinberg 1994 ; Evrard, Summers, fe Davis 1994| 



Anninos et al. 1994; Gnedin 1995). Therefore, one can ex- 



pect that consistency between those two different methods 
may not be automatically achieved. 

The cooling inconsistency was largely ignored in the 
pre vious numerical work except in a simulation by Katz et 
al. (1994), who have abandoned the cooling time condition 
in their three-dimensional simulation of galaxy formation for 
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all gas elements below 30, 000 K; this solution to the cooling 
inconsistency is acceptable as long as an equilibrium stan- 
dard cooling curve is used as a cooling function; in a sim- 
ulation which includes nonequilibrium time-dependent ion- 
ization, radiative transfer effects, cooling by heavy elements 
and/or by molecules, this approach would become inappro- 
priate, since a cooling function is then a function of position 
and time, and no characteristic temperature like 30, 000 K, 
below which the cooling time condition may be abandoned, 
can be specified a priori. 

This papers is composed as follows. In §2 I utilize the 
spherically symmetric Lagrangian hydrodynamic code to 
emphasize the importance and role of the cooling consis- 
tency condition, and in §3 I apply those results to cor- 
rect the existing cosmological hydrodynamic code based on 
the "Softened Lagrangian Hydrodynamics" method (SLH; 
Gnedin 1995) for the inconsistency between the hydrody- 
namic and cooling solvers. I conclude in §4 with brief dis- 



2 COOLING CATASTROPHE IN A 
SPHERICALLY SYMMETRIC CASE 

I first consider a collapse of a spherically symmetric self- 
gravitating gas cloud which is loosing its energy via radia- 
tive cooling. At the initial state the cloud is at rest, has a 
uniform temperature T = 1 (in dimensionless units) and the 
density profile: 



P = 



1 + r 2 



(1) 



where all quantities are dimensionless for simplicity, and the 
gravitational constant G = 1. The gas has a constant poly- 
tropic index 7 = 5/3, and equation of state is 

P = pT. 

The internal energy of the gas is lost due to radiative cooling 
at a rate 



(-) 

\ at J cool 



-pA(T), 



(2) 



where the cooling function A(T) has the following form: 

A(T) = 10" 2 T 5 exp(-T -4 ). (3) 

This cooling function allows for an analytical solution for the 
cooling evolution at a time interval short compared to the 
dynamical time and simultaneously mimics the steep cut-off 
at 10 4 K in the real cooling function. 

Apparently, the gas cloud would loose its energy due 
to radiative cooling and would collapse to the central sin- 
gularity in a finite time, the process known as the "cooling 
catastrophe". In order to simulate this process, the evolu- 
tion of the gas is followed numerically by a on e-dimensional 
Lagrangian code as described in Gnedin ( 199E| ) (adapted for 
the spherically symmetric case and combined with the ex- 
act gravity solver). Since the hydrodynamic solver is exactly 
Lagran gian, the gravitational consist ency condition becomes 
trivial (Gnedin & Bertschinger 1996). 



At each hydrodynamic time step t„ new values of the 
gas density p n , velocity v n , and temperature without cooling 
T n are computed using the values at the previous time step 
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Figure 1. The gas density [lower panel), velocity (middle panel), 
and temperature (upper panel) for the exact spherically symmet- 
ric collapse as a function of radius. Two moments are shown: the 
initial moment (thin line) and the final moment (bold line) de- 
fined as the moment when the central density has increased by 
12 orders of magnitude. The density profile is a perfect power-law 
with the index -2 (except near the center where numerical loss of 
accuracy occurs at very high densities). 



t n -i- Then, in order to account for radiative cooling, the 
following equation for the temperature is solved: 

^ = -(7-lKA(T) (4) 

in a time interval At — t n — t n -i with the initial condition 
T(t n -i) = T„0 This is the way the radiative cooling is 
incorporated in most of existing cosmological hydrodynamic 
codes. 

Figure [l] shows the gas density, velocity and tempera- 
ture as a function of radius at two different moments: the 
initial moment t = 0, shown as a thin solid line, and the mo- 
ment when the central density reached 10 12 , i.e. increased 
by 12 orders of the magnitude (the bold line). At this mo- 
ment the simulation with 256 zones encountered numerical 
problems at the center due to the round-off error in double 
precision calculations. The density distribution is a perfect 
power law with the index -2 across the 7 orders of magni- 
tude in radius and 14 orders of magnitude in density, and the 
gas temperature is almost constant, changing by less than 
0.2 dex (60%) across 7 orders of magnitude in radius. 

This results demonstrates both the success and failure 
of the purely Lagrangian gas dynamics. It definitely succeeds 
in achieving enormous dynamical range, and simultaneously 
it fails to carry on the solution beyond the collapse time. In- 
deed this is a correct result which tells one that the simple 
model one adopted is not sufficient to describe the behavior 



* With the cooling function (B 
analytically. 



this equation can be solved 
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Figure 2. The gas density (lower panel), velocity (middle panel), 
and temperature (upper panel) as a function of radius for the 
exact case (solid line) and the softened gravity case (dashed line). 
The softened case is shown for the time instant t = l.lt co i, where 
t co l is the collapse time of the exact solution. 

of a real physical system. Nevertheless, in a real three dimen- 
sional simulation, one would like to proceed further in time 
simulating other regions of the universe even if there has 
formed a real singularity in one point in space. Numerically 
this is achieved by softening the gravity solver, i.e. using a 
softened gravity law instead of the exact 1/r 2 Newtonian 
law. 

Let me now consider the cooling catastrophe for the 
softened gravity case. The standard Plummer softened grav- 
ity law is adopted with the softening length e = 0.01 and 
the same initial configuration is followed forward in time. 
In this case no real singularity forms, and the solution can 
be carried on in time indefinitely. Figure |^ shows the gas 
density, velocity and the temperature for the exact case as 
in FigJll with a solid line, and for the softened case with the 
dashed line. The softened solution is shown at a moment 
t = l.ltcoi, where t co i is the time of collapse for the exact 
solution. As one can expect, the velocity approaches zero, 
and the density and the temperature approach finite values 
at the center, with the core radius close to the softening 
length e. 

As one can expect, the softened solution provides a good 
approximation to the real solution on scales larger than sev- 
eral softening lengths. The temperature inside the softening 
length is only slightly lower than the exact solution and thus 
we may conclude that the softened solution is a good approx- 
imation to the exact solution. Here, however, the cooling 
consistency condition steps in. 

When the gravitational consistency condition is vio- 
lated, there is a simple basic law of physics - energy conser- 
vation - which is also violated, and it is then clear that the 
gravitational consistency condition has to assure the con- 
servation of energy. It is not quite so with the cooling con- 
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Figure 3. The gas cooling time (lower panel) and temperature 
(upper panel) for three cases: exact solution (solid line), softened 
gravity without cooling consistency condition (dashed line), and 
softened gravity with the cooling consistency condition (dotted 
line). The balance between heating and cooling is broken in the 
softened gravity solution which is manifested in the unphysically 
long cooling times; the balance is restored back in the consistent 
solution. 

sistency condition, since there is no basic physical law that 
is violated by the softened solution in Figj2| If, for exam- 
ple, one is satisfied with the mere existence of the gas globs 
described by the dashed line in Figji] in one's simulation, 
one can simply ignore any cooling inconsistency that may 
appear in the softened solution. However, if one does not re- 
strict himself to the mere existence of those gas globs, but is 
interested whether those objects actually cool and collapse 
to eventually form stars (as the exact solution does), the 
softened solution becomes unacceptable. In Figure ^ I show 
the cooling time and the temperature as a function of ra- 
dius for the exact solution (the solid line) and the softened 
solution (the dashed line; the dotted line is explained be- 
low) . Even if temperatures of two solutions differ by a small 
amount, the cooling function is so steep, that the cooling 
time of the softened solution exceed the cooling time of the 
exact solution by many orders of magnitude. At every point 
in the exact solution the cooling time tcooL is similar to the 
local dynamical time £dyn , 

icooL ~ £dyn- 

This condition is satisfied because the collapse proceeds at 
the dynamical time, and if the cooling time is much shorter 
than the dynamical time, the temperature will quickly de- 
crease, the cooling function will also decrease and the cool- 
ing time will increase until it is close to the dynamical time; 
if, inversely, the dynamical time is shorter than the cooling 
time, the temperature will increase adiabatically, the cool- 
ing function will increase and the cooling time will decrease 
until it is again approximately equal to the dynamical time. 
The collapse is therefore proceeding in the quasi-equilibrium 
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state, when adiabatic heating is approximately balanced by 
radiative cooling. 

The softened solution, however, grossly violates this 
quasi-equilibrium. Since the softened gravity prevents the 
collapse much beyond the resolution limit of the code, there 
is no adiabatic heating to balance cooling, the gas temper- 
ature will continue to decrease and the cooling time will 
continue to increase as long as the simulation progresses, so 
that at a sufficient late time moment t, 

tcooh ~ t ^> £dyn • 

Then, if one tries to estimate whether the object is actually 
cooling and collapsing by comparing its dynamical time to 
its cooling time, one would find that the cooling time is 
several orders of magnitude larger than the dynamical time, 
and would erroneously conclude that the object does not 
cool and, therefore, does not collapse! 

This conclusion stems from the fact that the hydrody- 
namic and cooling solvers are not consistent with each other 
in a sense that the cooling solver (^) knows nothing about 
whether the resolution limit of the code has been already 
reached or not. I can therefore immediately conclude that 
for a fully Lagrangian solver the cooling consistency condi- 
tion is trivial as is the gravitational consistency condition. 
However, in this case it is assumed that the fully Lagrangian 
solver obeys the 1/r 2 gravity law, whereas for the gravita- 
tional con sistency condition only poten tiality of the force is 
required (Gnedin & Bertschinger 1996). 



Let me now consider the softened case and address the 
question whether we can repair the deficiency of the cooling 
solver and make it consistent with the (self-gravitating) hy- 
drodynamic solver. Since the inconsistency appears as the 
"overcooling" problem when the resolution limit is reached, 
the obvious solution is to reduce the cooling rate at the res- 
olution limit of the code. I therefore replace the radiative 
cooling term (H) with the following expression: 



-Dp A(T), 



(5) 



(-) 

V dt J cool 

where D is the cooling consistency factor, < D < 1. Ob- 
viously, D should be a function of the softening length e, so 
that D = 1 when r 2> e and D — if r <C e. To demonstrate 
the effect of the cooling consistency condition, I adopt the 
following form for the cooling consistency factor D: 

D (.^)=(-^-s) , (6) 



where 



1, 

0.01, 
0.01(e/r 



if r < O.le, 
if r > e, and 
otherwise, 



(7) 



and 



if r < O.Ole, 
if r > O.le, and 
otherwise. 



(8) 



r 0.75, 
b = I 0.85, 

1 0.75 + 0.1(2 + log 10 (r/e)), 

The dotted line in Fig.^ shows the temperature and the cool- 
ing time for the softened case when the cooling consistency 
condition is incorporated. The cooling time agrees very well 
with the exact solution. The density and the velocity pro- 
files for the corrected case coincide with the uncorrected 
softened case almost precisely, so the cooling consistency 



condition does not change the density profile. Apparently, 
it is possible to achieve even better agreement by an ap- 
propriate choice of the cooling consistency factor D. The 
particular choice presented above serves merely to demon- 
strate the importance and effect of the cooling consistency 
condition. 



3 COOLING CONSISTENCY CONDITION IN 
THREE DIMENSIONS 

Let me first briefly discuss the rationale behind introducing 
the cooling consistency factor D. The equation for the tem- 
perature of a fluid element with a constant polytropic index 
7 can be written as: 



dF = (7 - 1} 



(9) 



where d/dt denotes the Lagrangian time-derivative, and the 
cooling function A is not required to be a function of T 
only. For a cooling catastrophe, the dynamical time is close 
to the cooling time, which implies that the two terms in 
square brackets almost cancel each other and the collapse is 
almost isothermal: 



(10) 



When equation is implemented numerically, the time 
derivative becomes numerical time derivative for a finite res- 
olution numerical scheme, which I emphasize by adding a 
subscript N to it: 



djyT 
dt 



(7-1) 



Td N p 
p dt 



pA(T,p) 



(11) 



The main difference between equations Q and ( pd| ) is that 
at the resolution limit of the code the density stops changing, 



d N p 



(res. limit) 



and equation 

djyT 
dt 



►0, 

becomes: 



= -( 7 -l)pA(T,p), 



(12) 



(13) 



implying that the cooling time will increase in proportion to 
the physical time t and not the dynamical time tDYN- The 
equation ([nj) can be "repaired" by introducing the cooling 
consistency factor D, 



djyT 
dt 



(7-1) 



Tdjvp 
p dt 



D P K{T,p) 



Now, by requiring that 
djvp/df 



D = 



dp/dt 



(14) 



(15) 



I restore the quasi-equilibrium character of the collapse (eq. 
flTo||). The only disadvantage of this approach is that the 
exact cooling consistency condition ( p^[ ) can never be sat- 
isfied since the exact time derivate of the density, dp/dt is 
unknown, and more that that, it is what one tries to approx- 
imate with the numerical derivative djvp/dt. The cooling 
consistency condition can therefore be satisfied only approx- 
imately. 
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Figure 4. The gas cooling time vs dynamical time (in arbitrary 
units) for realistic three-dimensional cosmological simulations of a 
standard CDM model: a 160 dynamical range simulation with 5% 
baryons and without the cooling consistency condition (upper left 
panel); a 160 dynamical range simulation with 50% baryons and 
without the cooling consistency condition (lower left panel); a 160 
dynamical range simulation with 5% baryons and with the cooling 
consistency condition (upper right panel); a 640 dynamical range 
simulation with 5% baryons and without the cooling consistency 
condition (lower right panel); Again, the balance between heating 
and cooling is broken without the consistency condition indepen- 
dently of the baryonic fraction; the balance is restored by the 
cooling consistency condition or by increasing the resolution of a 
simulation. In the latter case however, cooling and heating rates 
get out of balance again at the (finer) resolution limit of a higher 
resolution simulation. 



Reducing cooling is not the only way to achieve cooling 
consistency. It is possible, for example, to introduce addi- 
tional heating terms in equation ( pd| ) to account for the loss 
of adiabatic heating near the resolution limit of the code. 
Therefore, there exist many different ways to satisfy the 
cooling consistency condition, and the cooling consistency 
factor D is only one of them. 

Let me now turn to a three-dimensional ca se. I use the 
SLH-P 3 M code (|Gnedin fc Bertschinger 1996|) to demon- 



strate the effect of the cooling consistency condition on the 
thermal evolution of the cosmic gas. I adopt the standard 
CDM model as a testbed for my demonstration. The cosmo- 
logical parameters are fixed so that f2o = 1-0, h = 0.5, and 
I fix the i nitial power spectrum by using the BBKS transfer 
function ( Bardeen et al. 1986 ) and adopting as = 1.0. 

First, two 32 3 runs with the softening parameter of 1/5 
(the dynamical range of 160) and without any consistency 
condition were performed, the first one using Qb = 0.05, 
and, therefore, having 5% baryons, and the second one us- 



ing fib = 0.5 (50% baryons) with otherwise identical initial 
conditions. Figure ^ shows cooling and dynamical times for 
all fluid elements with overdensities in excess of 100 and 
temperatures between 10 3 ' 5 K and 10 4 ' 5 K at z = 6 in a 
2/i _1 Mpc box (at this moment most of the gas has not yet 
been shocked to temperatures significantly exceeding 10 4 K) 
for those two runs in the upper left and lower left panels re- 
spectively. It is apparent that the cooling time in the highest 
density regions is very long, exceeding the dynamical time by 
up to two orders of magnitude, and the balance between the 
cooling and heating rates is severely broken - the manifesta- 
tion of the cooling inconsistency. Were those fluid elements 
tested on the condition of collapse, one would find that they 
cannot cool, and, therefore, cannot collapse, which is appar- 
ently nonphysical conclusion. As one can expect, the cool- 
ing inconsistency does not depend on the baryonic fraction 
((and, hence, on the dark matter fraction) since equation 
( pT| ) contains no reference to the dark matter fraction. 

I can now introduce the SLH cooling consistency con- 
dition by defining the cooling consistency factor D in the 
following way: 



D = 1 - min(cri, 0-2,0-3), 



(16) 



where Oj are eigenvalues of the deformation tensor o 1J which 
defines the resolution limit of the SLH c ode (for exact defi- 
nitions see Gnedin & Bertschinger ( 1996] )). The cooling con- 
sistency condition in this form gradually reduces the total 
cooling rate as the gas element followed by the SLH code 
approaches the code resolution limit. 

Since the SLH code is not exactly Lagrangian, and, in 
general, the flow is not spherically symmetric, an expres- 
sion for D found in the spherically symmetric case (eq. Q) 
cannot be applied in a three-dimensional case exactly. How- 
ever, in the spherically symmetric case, the expression ( |16 ) 
reduces to equation (g) with a — b = 1. While corrections (17 ) 
and (^|) improve the accuracy of the approximate solution, 
they cannot be easily generalized for a three-dimensional 
case. The three-dimensional consistency condition ([l6]) is, 
therefore, less accurate that the three-dimensional approx- 
imation given by (jfjj); further improving upon the accuracy 
of the SLH consistency condition is beyond the frame of this 
paper. 

I can now test the SLH consistency condition (jlq). The 
distribution of cooling and dynamical times for a 32^ simu- 
lation with the 1/5 softening and Q(, = 0.05 with the cooling 
consistency condition as given by equation jl6[ ) is shown in 
Fig.^ at the upper right panel. One can note now that in 
the run with the consistency condition, the cooling times 
are in a reasonable (yet not perfect, since the cooling con- 
sistency condition ( |l6[ ) is only an approximate one) agree- 
ment with the dynamical times for all fluid elements with 
high density. For those fluid elements the cooling time is of 
the order of their dynamical time and they are apparently 
cooling and collapsing. One must remember, of course, that 
in the simulation most of those fluid elements have already 
reached their resolution limit and their density does not sig- 
nificantly increase with time, but the relationship between 
their cooling and dynamical times mimics that of truly cool- 
ing and collapsing gas clouds. A simulation with the cooling 
consistency condition and Qb = 0.5 (50% baryons) produces 
nearly identical distribution of cooling and dynamical times 
and is not shown here due to the space limitations. 



© 0000 RAS, MNRAS 000, 000-000 



6 Gnedin 



Finally, the lower right panel of Fig.^ shows the distri- 
bution of cooling and dynamical times for a 64 3 simulation 
with the softening parameter of 1/10 (the dynamical range 
of 640) with the same 2/i _1 Mpc box and Qt — 0.05 with- 
out the cooling consistency condition (for clarity only every 
eights fluid element is shown) . Since this simulation has four 
times higher spatial resolution than the simulation shown in 
the upper left panel of Fig.^, the balance between cooling 
and heating extends to about a factor of 64 higher densities, 
or, equivalently, to about a factor of 8 smaller dynamical 
times. One can easily see that the simulation with the consis- 
tency condition (the upper right panel) mimics very closely 
the high resolution simulation. At the highest densities the 
balance is again broken, as manifested by the cooling time 
becoming larger than the dynamical time as the dynamical 
time decreases. The effect is not as dramatic as in the up- 
per left panel, and the gas densities in the high resolution 
simulation are only slightly larger than in the low resolution 
simulations because the final epoch (z — 6) was chosen such 
that in the high resolution simulation only a small number 
of all fluid elements have actually reached their resolution 
limit; this makes the high resolution simulation closely mim- 
icking an imaginary exact (infinite resolution) case. 
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4 CONCLUSIONS 

I have shown on simple spherically symmetric and realistic 
three-dimensional examples how the inconsistency between 
the finite resolution self-gravitating hydrodynamic solver 
and the radiative cooling term may arise. For a efficiently 
cooling and collapsing gas clouds the cooling time approxi- 
mately balances the dynamical time so that the collapse is 
occurring in the quasi-equilibrium between adiabatic heat- 
ing and radiative cooling. 

In simulations, the finite resolution leads to decrease in 
adiabatic heating when a fluid elements approaches the res- 
olution limit of a numerical code, and the balance between 
heating and coolings breaks down. This leads to an "over- 
cooling" problem, when the cooling time in the fluid element 
becomes significantly longer than the dynamical time, and 
the fluid element appears cooling very inefficiently. 

The contradiction is eliminated when the cooling con- 
sistency condition is introduced, which reduces the cooling 
rate in proportion to the reduction in the heating rate. I 
present the exact cooling consistency condition, which how- 
ever cannot be realized in practice since it depends on the 
unknown true solution which a simulations tries to approxi- 
mate. However, approximate cooling consistency conditions 
can be introduced both in a spherically symmetric and re- 
alistic three-dimensional cases which are "good enough" in 
a sense that the resultant approximate numerical solution 
mimics the exact solution in balancing the heating and cool- 
ing rates. 
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